#*******************************************************************************************
# Plots for the manuscript:
# Germany’s decision to phase out coal by 2038 lags behind citizens‘ timing preferences
# Author: Adrian Rinscheid
#*******************************************************************************************

setwd("/Users/Someone/Dropbox")

rm(list=ls())
library("ggplot2")

# Figure 2 | Average effects of policy attributes on respondents’ preference for a coal phase-out.

d <- read.table("coal_brating40.txt",comment.char="")
ffilename <- "Fig.2"

# group vars
Enddate     <- paste(c(1:2,"3b",4),".attime",sep="")
Costs     <- paste(c("1b",2:4),".atcost",sep="")
Lostjobs    <- paste(c("1b",2:4),".atlostjobs",sep="")
Newjobs   <- paste(c("1b",2:4),".atnewjobs",sep="")
Measures       <- paste(c("1b",2:5),".atmeasures",sep="")


CIs <- function(d){
  colnames(d)[1:2] <- c("pe","se")
  d$upper <-d$pe + 1.96*d$se
  d$lower <-d$pe - 1.96*d$se
  return(d)
}
d<- CIs(d)
d$var <- rownames(d)

FillGroup <- function(d){
  
  d$gruppe <- NA
  d$gruppe[d$var %in% Enddate]            <- "End date"
  d$gruppe[d$var %in% Costs]              <- "Annual costs"
  d$gruppe[d$var %in% Lostjobs]           <- "Lost jobs"
  d$gruppe[d$var %in% Newjobs]            <- "New jobs"
  d$gruppe[d$var %in% Measures]           <- "Supportive measures"
  
  return(d)
}

d <- FillGroup(d)
d$order <- 1:nrow(d)

GetLabels <- function(d){
  offset <- c("   ")
  
  d$var[d$var %in% Enddate] <- paste(offset,c("2025","2030",
                                              "2040","2100"))
  
  d$var[d$var %in% Costs] <- paste(offset,c("EUR 0", "EUR 6", "EUR 12",
                                            "EUR 18"))
  
  d$var[d$var %in% Lostjobs] <- paste(offset,c("-5,000","-10,000","-15,000",
                                               "-20,000"))
  
  d$var[d$var %in% Newjobs] <- paste(offset,c("5,000","10,000","15,000",
                                              "20,000"))
  
  d$var[d$var %in% Measures] <- paste(offset,c("Expansion of renewables",
                                               "Support for new businesses",
                                               "Investment in modern infrastructure",
                                               "Research & development",
                                               "Training & early retirement"))
  return(d)
}

d <- GetLabels(d)

# bring in sublabels           
d$order <- 1:nrow(d)
dd <- data.frame(var= c("End date:",
                        " ",
                        "Annual costs:",
                        "  ",
                        "Lost jobs:",
                        "   ",
                        "New jobs:",
                        "    ",
                        "Supportive measures:"
),order=c(.5,4.1,4.5,8.1,8.5,12.1,12.5,16.1,16.5),
pe=1,se=1,upper=1,lower=1,gruppe=NA)
d <- rbind(d,dd)
d <-d[order(d$order),]
d$var <- factor(d$var,levels=unique(d$var)[length(d$var):1])

yylab  <- c("Phase-out Support (binary indicator)")

p = ggplot(d,aes(y=pe,x=var))#,colour=gruppe))
p = p + coord_flip(ylim = c(-.3, .3))  
p = p + geom_hline(yintercept = 0,size=.5,colour="black",linetype="dotted") 
p = p +geom_pointrange(aes(ymin=lower,ymax=upper,width=.4),position="dodge",size=.6)
p = p + scale_y_continuous(name=yylab,breaks=round(seq(-.4,.4,.2),1),labels=c("-.4","-.2","0",".2",".4"))
p = p + scale_colour_discrete("Attribute:") + scale_x_discrete(name="")
print(p)

theme_bw1 <- function(base_size = 13, base_family = "") {
  theme_grey(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.text.x =       element_text(size = base_size, colour = "black",  hjust = .5 , vjust=1),
      axis.text.y =       element_text(size = base_size , colour = "black", hjust = 0 , vjust=.5 ), # changes position of X axis text
      axis.ticks =        element_line(colour = "grey50"),
      axis.title.y =      element_text(size = base_size,angle=90,vjust=.01,hjust=.1),
      legend.position = "none"
    )
}

dev.off()
pdf(paste(ffilename,".pdf",sep=""),width=10,height=10) # for paper
p = p  + theme_bw1()
print(p)
dev.off()


# Figure 3 | Average effects of timing attribute on respondents’ preference for a coal-phase out by party identification and perceived scientific consensus.

# Fig. 3a: party identification

Sys.setlocale("LC_ALL", 'en_US.UTF-8')#allows for ä, ö, ü in the plot

library(ggplot2)
library(foreign)

# Read in data
d1 <- read.dta("coal_parties440.dta")
d1 <- d1[d1[,1]!="_cons",]
d1 <- d1[1:4,]
d1$group <- rownames(d1)
d1$model <- "CDU"
d1$model_id <- "1"

#Create Vectors for coefs and standard errors for each model, and variable names
coef.vec.1<-d1$estimate
se.vec.1 <- d1$stderr

d2 <- read.dta("coal_parties240.dta")
d2 <- d2[d2[,1]!="_cons",]
d2 <- d2[1:4,]
d2$group <- rownames(d2)
d2$model <- "SPD"
d2$model_id <- "2"

d3 <- read.dta("coal_parties340.dta")
d3 <- d3[d3[,1]!="_cons",]
d3 <- d3[1:4,]
d3$group <- rownames(d3)
d3$model <- "Die Linke"
d3$model_id <- "3"

d4 <- read.dta("coal_parties140.dta")
d4 <- d4[d4[,1]!="_cons",]
d4 <- d4[1:4,]
d4$group <- rownames(d4)
d4$model <- "Die Grünen"
d4$model_id <- "4"

d5 <- read.dta("coal_parties540.dta")
d5 <- d5[d5[,1]!="_cons",]
d5 <- d5[1:4,]
d5$group <- rownames(d5)
d5$model <- "CSU"
d5$model_id <- "5"

d6 <- read.dta("coal_parties640.dta")
d6 <- d6[d6[,1]!="_cons",]
d6 <- d6[1:4,]
d6$group <- rownames(d6)
d6$model <- "FDP"
d6$model_id <- "6"

d7 <- read.dta("coal_parties740.dta")
d7 <- d7[d7[,1]!="_cons",]
d7 <- d7[1:4,]
d7$group <- rownames(d7)
d7$model <- "AFD"
d7$model_id <- "7"


#Create Vectors for coefs and standard errors for each model, and variable names
coef.vec.2<-d2$estimate
se.vec.2 <- d2$stderr
coef.vec.3<-d3$estimate
se.vec.3 <- d3$stderr
coef.vec.4<-d4$estimate
se.vec.4 <- d4$stderr
coef.vec.5<-d5$estimate
se.vec.5 <- d5$stderr
coef.vec.6<-d6$estimate
se.vec.6 <- d6$stderr
coef.vec.7<-d7$estimate
se.vec.7 <- d7$stderr

var.names <- c("2025", "2030", "2040", "2100"
)

y.axis <- length(var.names): 1 
adjust <- .04 

pdf("Fig.3a.pdf", encoding = "ISOLatin9.enc",width=9,height=7.5)

layout(matrix(c(2,1),1,2), 
       widths = c(.2, 5))

par(mar=c(4,4.5,.5,.2), lheight = .8)

plot(coef.vec.1, y.axis, type = "p", axes = F, xlab="", ylab = "", pch = 17, cex = 1.7, 
     xlim = c(min((coef.vec.2-qnorm(.975)*se.vec.2 -.1), (coef.vec.3-qnorm(.975)*se.vec.3 -.1), (coef.vec.4-qnorm(.975)*se.vec.4 -.1), (coef.vec.5-qnorm(.975)*se.vec.5 -.1), (coef.vec.6-qnorm(.975)*se.vec.6 -.1), (coef.vec.7-qnorm(.975)*se.vec.7 -.1), na.rm = T), #set xlims at mins and maximums (from both models) of confidence intervals, plus .1 to leave room at ends of plots
              max((coef.vec.2+qnorm(.975)*se.vec.2 +.1), (coef.vec.3+qnorm(.975)*se.vec.3 +.1), (coef.vec.4+qnorm(.975)*se.vec.4 +.1), (coef.vec.5+qnorm(.975)*se.vec.5 +.1), (coef.vec.6+qnorm(.975)*se.vec.6 +.1), (coef.vec.7+qnorm(.975)*se.vec.7 +.1), na.rm = T)),  
     ylim = c(min(y.axis), max(y.axis)), main = "")

title(xlab="Phase-out Support (binary indicator)", cex.lab = 1.5)
title(ylab="End date", cex.lab = 1.5)
axis(1, at =seq(-.6,.3,.1), label = seq(-.6,.3,.1), cex.axis = 1.5)
axis(2, at = y.axis, label = var.names, las = 1, tick = T, cex.axis =1.5)
abline(h = y.axis, lty = 2, lwd = .5, col = "grey")

segments(coef.vec.1-qnorm(.975)*se.vec.1, y.axis, coef.vec.1+qnorm(.975)*se.vec.1, y.axis, lwd =  1.3)
abline(v=0, lty = 2) 

segments(coef.vec.2-qnorm(.975)*se.vec.2, y.axis+adjust*1.2, coef.vec.2+qnorm(.975)*se.vec.2, y.axis+adjust*1.2, col="violet", lwd =  1.3)
points(coef.vec.2, y.axis+adjust*1.2, pch = 17, col="violet", cex = 1.7, bg = "white" ) 

segments(coef.vec.3-qnorm(.975)*se.vec.3, y.axis+adjust*.6, coef.vec.3+qnorm(.975)*se.vec.3, y.axis+adjust*.6, col="purple", lwd =  1.3)
points(coef.vec.3, y.axis+adjust*.6, pch = 19, col="purple", cex = 1.7, bg = "white" ) 

segments(coef.vec.4-qnorm(.975)*se.vec.4, y.axis+adjust*1.8, coef.vec.4+qnorm(.975)*se.vec.4, y.axis+adjust*1.8, col="darkgreen", lwd =  1.3)
points(coef.vec.4, y.axis+adjust*1.8, pch = 18, col="darkgreen", cex = 1.9, bg = "white" ) 

segments(coef.vec.5-qnorm(.975)*se.vec.5, y.axis-adjust*.6, coef.vec.5+qnorm(.975)*se.vec.5, y.axis-adjust*.6, col="gray", lwd =  1.3)
points(coef.vec.5, y.axis-adjust*.6, pch = 15, col="gray", cex = 1.7, bg = "white" ) 

segments(coef.vec.6-qnorm(.975)*se.vec.6, y.axis-adjust*1.2, coef.vec.6+qnorm(.975)*se.vec.6, y.axis-adjust*1.2, col="orange", lwd =  1.3)
points(coef.vec.6, y.axis-adjust*1.2, pch = 18, col="orange", cex = 1.9, bg = "white" ) 

segments(coef.vec.7-qnorm(.975)*se.vec.7, y.axis-adjust*1.8, coef.vec.7+qnorm(.975)*se.vec.7, y.axis-adjust*1.8, col="darkblue", lwd =  1.3)
points(coef.vec.7, y.axis-adjust*1.8, pch = 16, col="darkblue", cex = 1.7, bg = "white" ) 

text(-.36, 3.9, "Die Grünen (n=138)", adj = 0,cex  = 1.3)
points(-.39, 3.9, pch = 18, col="darkgreen", cex  = 1.9)
text(-.36, 3.76, "SPD (n=325)", adj = 0,cex  = 1.3)
points(-.39, 3.76, pch = 17, col="violet", cex  = 1.7)
text(-.36, 3.62, "Die Linke (n=160)", adj = 0,cex  = 1.3)
points(-.39, 3.62, pch = 19, col="purple", cex  = 1.7)
text(-.36, 3.48, "CDU (n=321)", adj = 0,cex  = 1.3)
points(-.39, 3.48, pch = 17, cex  = 1.7)
text(-.36, 3.34, "CSU (n=85)", adj = 0,cex  = 1.3)
points(-.39, 3.34, pch = 15, col="gray", cex  = 1.7)
text(-.36, 3.2, "FDP (n=81)", adj = 0,cex  = 1.3)
points(-.39, 3.2, pch = 18, col="orange", cex  = 1.9)
text(-.36, 3.06, "AFD (n=129)", adj = 0,cex  = 1.3)
points(-.39, 3.06, pch = 16, col="darkblue", cex  = 1.7)
dev.off()


# Fig. 3b: perceived scientific consensus.

Sys.setlocale("LC_ALL", 'en_US.UTF-8')#allows for ä, ö, ü in the plot

library(ggplot2)
library(foreign)

# Read in data
d1 <- read.dta("coal_consensus240.dta")
d1 <- d1[d1[,1]!="_cons",]
d1 <- d1[1:4,]
d1$group <- rownames(d1)
d1$model <- "50-69%"
d1$model_id <- "1"

#Create Vectors for coefs and standard errors for each model, and variable names
coef.vec.1<-d1$estimate
se.vec.1 <- d1$stderr

d2 <- read.dta("coal_consensus140.dta")
d2 <- d2[d2[,1]!="_cons",]
d2 <- d2[1:4,]
d2$group <- rownames(d2)
d2$model <- "0-49%"
d2$model_id <- "2"

d3 <- read.dta("coal_consensus340b.dta")
d3 <- d3[d3[,1]!="_cons",]
d3 <- d3[1:4,]
d3$group <- rownames(d3)
d3$model <- "70-84%"
d3$model_id <- "3"

d4 <- read.dta("coal_consensus440b.dta")
d4 <- d4[d4[,1]!="_cons",]
d4 <- d4[1:4,]
d4$group <- rownames(d4)
d4$model <- "85-100%"
d4$model_id <- "4"


#Create Vectors for coefs and standard errors for each model, and variable names
coef.vec.2<-d2$estimate
se.vec.2 <- d2$stderr
coef.vec.3<-d3$estimate
se.vec.3 <- d3$stderr
coef.vec.4<-d4$estimate
se.vec.4 <- d4$stderr

var.names <- c("2025", "2030", "2040", "2100"
)

y.axis <- length(var.names): 1 
adjust <- .04 

pdf("Fig.3b.pdf", encoding = "ISOLatin9.enc",width=9,height=7.5)

layout(matrix(c(2,1),1,2), 
       widths = c(.2, 5))

par(mar=c(4,4.5,.5,.2), lheight = .8)

plot(coef.vec.1, y.axis, type = "p", axes = F, xlab="", ylab = "", pch = 18, col="orange", cex = 1.9, 
     xlim = c(min((coef.vec.2-qnorm(.975)*se.vec.2 -.1), (coef.vec.3-qnorm(.975)*se.vec.3 -.1), (coef.vec.4-qnorm(.975)*se.vec.4 -.1), na.rm = T), 
              max((coef.vec.2+qnorm(.975)*se.vec.2 +.1), (coef.vec.3+qnorm(.975)*se.vec.3 +.1), (coef.vec.4+qnorm(.975)*se.vec.4 +.1), na.rm = T)),  
     ylim = c(min(y.axis), max(y.axis)), main = "")

title(xlab="Phase-out Support (binary indicator)", cex.lab = 1.5)
title(ylab="End date", cex.lab = 1.5)
axis(1, at =seq(-.6,.3,.1), label = seq(-.6,.3,.1), cex.axis = 1.5)
axis(2, at = y.axis, label = var.names, las = 1, tick = T, cex.axis =1.5)
abline(h = y.axis, lty = 2, lwd = .5, col = "grey")

segments(coef.vec.1-qnorm(.975)*se.vec.1, y.axis, coef.vec.1+qnorm(.975)*se.vec.1, y.axis, col="orange", lwd =  1.3)
abline(v=0, lty = 2) 

segments(coef.vec.2-qnorm(.975)*se.vec.2, y.axis+adjust*.85, coef.vec.2+qnorm(.975)*se.vec.2, y.axis+adjust*.85, col="red", lwd =  1.3)
points(coef.vec.2, y.axis+adjust*.85, pch = 17, col="red", cex = 1.7, bg = "white" ) 

segments(coef.vec.3-qnorm(.975)*se.vec.3, y.axis-adjust*.85, coef.vec.3+qnorm(.975)*se.vec.3, y.axis-adjust*.85, col="lightblue", lwd =  1.3)
points(coef.vec.3, y.axis-adjust*.85, pch = 15, col="lightblue", cex = 1.7, bg = "white" ) 

segments(coef.vec.4-qnorm(.975)*se.vec.4, y.axis-adjust*1.7, coef.vec.4+qnorm(.975)*se.vec.4, y.axis-adjust*1.7, col="darkblue", lwd =  1.3)
points(coef.vec.4, y.axis-adjust*1.7, pch = 16, col="darkblue", cex = 1.8, bg = "white" ) 


text(-.32, 3.9, "00-49% perceived c. (n=368)", adj = 0,cex  = 1.3)
points(-.355, 3.9, pch = 17, col="red", cex  = 1.7)
text(-.32, 3.76, "50-69% perceived c. (n=513)", adj = 0,cex  = 1.3)
points(-.355, 3.76, pch = 18, col="orange", cex  = 1.9)
text(-.32, 3.62, "70-89% perceived c. (n=768)", adj = 0,cex  = 1.3)
points(-.355, 3.62, pch = 15, col="lightblue", cex  = 1.7)
text(-.32, 3.48, "90-100% perceived c. (n=390)", adj = 0,cex  = 1.3)
points(-.355, 3.48, pch = 16, col="darkblue", cex  = 1.8)
dev.off()


# Figure 4 | Average effects of timing attribute on respondents’ preference for a coal-phase out in Rhineland and Lusatia.

# Fig. 4a: Rhineland

Sys.setlocale("LC_ALL", 'en_US.UTF-8')#allows for ä, ö, ü in the plot

library(ggplot2)
library(foreign)

# Read in data
d1 <- read.dta("coal_brating240.dta")
d1 <- d1[d1[,1]!="_cons",]
d1 <- d1[1:4,]
d1$group <- rownames(d1)
d1$model <- "Rhineland_total"
d1$model_id <- "1"

#Create Vectors for coefs and standard errors for each model, and variable names
coef.vec.1<-d1$estimate
se.vec.1 <- d1$stderr

d2 <- read.dta("coal_worksknows240.dta")
d2 <- d2[d2[,1]!="_cons",]
d2 <- d2[1:4,]
d2$group <- rownames(d2)
d2$model <- "Rhineland_knows"
d2$model_id <- "2"

d3 <- read.dta("coal_notworksknows240.dta")
d3 <- d3[d3[,1]!="_cons",]
d3 <- d3[1:4,]
d3$group <- rownames(d3)
d3$model <- "Rhineland_notknows"
d3$model_id <- "3"


#Create Vectors for coefs and standard errors for each model, and variable names
coef.vec.2<-d2$estimate
se.vec.2 <- d2$stderr
coef.vec.3<-d3$estimate
se.vec.3 <- d3$stderr

var.names <- c("2025", "2030", "2040", "2100"
)


y.axis <- length(var.names): 1 
adjust <- .02 

pdf("Fig.4a.pdf", encoding = "ISOLatin9.enc",width=9,height=7.5)

layout(matrix(c(2,1),1,2), 
       widths = c(.2, 5))

par(mar=c(4,4.5,.5,.2), lheight = .8)

plot(coef.vec.1, y.axis+adjust*3, type = "p", axes = F, xlab="", ylab = "", pch = 17, cex = 1.7, 
     xlim = c(min((coef.vec.2-qnorm(.975)*se.vec.2 -.1), (coef.vec.3-qnorm(.975)*se.vec.3 -.1), na.rm = T), 
              max((coef.vec.2+qnorm(.975)*se.vec.2 +.1), (coef.vec.3+qnorm(.975)*se.vec.3 +.1), na.rm = T)),  
     ylim = c(min(y.axis), max(y.axis)), main = "")

title(xlab="Phase-out Support (binary indicator)", cex.lab = 1.5)
title(ylab="End date", cex.lab = 1.5)
axis(1, at =seq(-.4,.3,.1), label = seq(-.4,.3,.1), cex.axis = 1.5)
axis(2, at = y.axis, label = var.names, las = 1, tick = T, cex.axis =1.5)
abline(h = y.axis, lty = 2, lwd = .5, col = "grey")

segments(coef.vec.1-qnorm(.975)*se.vec.1, y.axis+adjust*3, coef.vec.1+qnorm(.975)*se.vec.1, y.axis+adjust*3, lwd =  1.3)
abline(v=0, lty = 2) 

segments(coef.vec.2-qnorm(.975)*se.vec.2, y.axis-adjust*1, coef.vec.2+qnorm(.975)*se.vec.2, y.axis-adjust*1, col="red", lwd =  1.3)
points(coef.vec.2, y.axis-adjust*1, pch = 18, col="red", cex = 2, bg = "white" ) 

segments(coef.vec.3-qnorm(.975)*se.vec.3, y.axis-adjust*4, coef.vec.3+qnorm(.975)*se.vec.3, y.axis-adjust*4, col="darkblue", lwd =  1.3)
points(coef.vec.3, y.axis-adjust*4, pch = 15, col="darkblue", cex = 1.7, bg = "white" ) 


text(-.255, 3.9, "Rhineland, total (n=491)", adj = 0,cex=1.3)
points(-.27, 3.9, pch = 17, cex  = 1.7)
text(-.255, 3.78, "working in coal industry", adj = 0,cex=1.3)
points(-.27, 3.78, pch = 18, col="red", cex=2)
text(-.255, 3.66, "   or having acquaintances (n=191)", adj = 0,cex=1.3)
points(-.27, 3.66, pch = 18, col="white", cex=.1)
text(-.255, 3.54, "not working in coal industry,", adj = 0,cex=1.3)
points(-.27, 3.54, pch = 15, col="darkblue", cex=1.7)
text(-.255, 3.42, "   no acquaintances (n=300)", adj = 0,cex=1.3)
points(-.27, 3.42, pch = 15, col="white", cex=.1)
dev.off()



# Fig. 4b: Lusatia

Sys.setlocale("LC_ALL", 'en_US.UTF-8')#allows for ä, ö, ü in the plot

library(ggplot2)
library(foreign)

# Read in data
d1 <- read.dta("coal_brating340.dta")
d1 <- d1[d1[,1]!="_cons",]
d1 <- d1[1:4,]
d1$group <- rownames(d1)
d1$model <- "Lusatia_total"
d1$model_id <- "1"

#Create Vectors for coefs and standard errors for each model, and variable names
coef.vec.1<-d1$estimate
se.vec.1 <- d1$stderr

d2 <- read.dta("coal_worksknows340.dta")
d2 <- d2[d2[,1]!="_cons",]
d2 <- d2[1:4,]
d2$group <- rownames(d2)
d2$model <- "Lusatia_knows"
d2$model_id <- "2"

d3 <- read.dta("coal_notworksknows340.dta")
d3 <- d3[d3[,1]!="_cons",]
d3 <- d3[1:4,]
d3$group <- rownames(d3)
d3$model <- "Lusatia_notknows"
d3$model_id <- "3"


#Create Vectors for coefs and standard errors for each model, and variable names
coef.vec.2<-d2$estimate
se.vec.2 <- d2$stderr
coef.vec.3<-d3$estimate
se.vec.3 <- d3$stderr

var.names <- c("2025", "2030", "2040", "2100"
)


y.axis <- length(var.names): 1 
adjust <- .02 

pdf("Fig.4b.pdf", encoding = "ISOLatin9.enc",width=9,height=7.5)

layout(matrix(c(2,1),1,2), 
       widths = c(.2, 5))

par(mar=c(4,4.5,.5,.2), lheight = .8)

plot(coef.vec.1, y.axis+adjust*3, type = "p", axes = F, xlab="", ylab = "", pch = 17, cex = 1.7, 
     xlim = c(min((coef.vec.2-qnorm(.975)*se.vec.2 -.1), (coef.vec.3-qnorm(.975)*se.vec.3 -.1), na.rm = T),
              max((coef.vec.2+qnorm(.975)*se.vec.2 +.1), (coef.vec.3+qnorm(.975)*se.vec.3 +.1), na.rm = T)), 
     ylim = c(min(y.axis), max(y.axis)), main = "")

title(xlab="Phase-out Support (binary indicator)", cex.lab = 1.5)
title(ylab="End date", cex.lab = 1.5)
axis(1, at =seq(-.4,.3,.1), label = seq(-.4,.3,.1), cex.axis = 1.5)
axis(2, at = y.axis, label = var.names, las = 1, tick = T, cex.axis =1.5)
abline(h = y.axis, lty = 2, lwd = .5, col = "grey")

segments(coef.vec.1-qnorm(.975)*se.vec.1, y.axis+adjust*3, coef.vec.1+qnorm(.975)*se.vec.1, y.axis+adjust*3, lwd =  1.3)
abline(v=0, lty = 2) 

segments(coef.vec.2-qnorm(.975)*se.vec.2, y.axis-adjust*1, coef.vec.2+qnorm(.975)*se.vec.2, y.axis-adjust*1, col="red", lwd =  1.3)
points(coef.vec.2, y.axis-adjust*1, pch = 18, col="red", cex = 2, bg = "white" ) 

segments(coef.vec.3-qnorm(.975)*se.vec.3, y.axis-adjust*4, coef.vec.3+qnorm(.975)*se.vec.3, y.axis-adjust*4, col="darkblue", lwd =  1.3)
points(coef.vec.3, y.axis-adjust*4, pch = 15, col="darkblue", cex = 1.7, bg = "white" ) 


text(-.21, 3.9, "Lusatia, total (n=473)", adj = 0,cex=1.3)
points(-.225, 3.9, pch = 17, cex  = 1.7)
text(-.21, 3.78, "working in coal industry", adj = 0,cex=1.3)
points(-.225, 3.78, pch = 18, col="red", cex=2)
text(-.21, 3.66, "   or having acquaintances (n=258)", adj = 0,cex=1.3)
points(-.225, 3.66, pch = 18, col="white", cex=.1)
text(-.21, 3.54, "not working in coal industry, ", adj = 0,cex=1.3)
points(-.225, 3.54, pch = 15, col="darkblue", cex=1.7)
text(-.21, 3.42, "   no acquaintances (n=215)", adj = 0,cex=1.3)
points(-.225, 3.42, pch = 15, col="white", cex=.1)
dev.off()


# Supplementary Figure 2 | Average effects of policy attributes on respondents’ preference for a coal phase-out, based on the forced choice outcome.

rm(list=ls())
library("ggplot2")

d <- read.table("coal_selected40.txt",comment.char="")
ffilename <- "Fig.S2"

# group vars
Enddate     <- paste(c(1:2,"3b",4),".attime",sep="")
Costs     <- paste(c("1b",2:4),".atcost",sep="")
Lostjobs    <- paste(c("1b",2:4),".atlostjobs",sep="")
Newjobs   <- paste(c("1b",2:4),".atnewjobs",sep="")
Measures       <- paste(c("1b",2:5),".atmeasures",sep="")


CIs <- function(d){
  colnames(d)[1:2] <- c("pe","se")
  d$upper <-d$pe + 1.96*d$se
  d$lower <-d$pe - 1.96*d$se
  return(d)
}
d<- CIs(d)
d$var <- rownames(d)

FillGroup <- function(d){
  
  d$gruppe <- NA
  d$gruppe[d$var %in% Enddate]            <- "End date"
  d$gruppe[d$var %in% Costs]              <- "Annual costs"
  d$gruppe[d$var %in% Lostjobs]           <- "Lost jobs"
  d$gruppe[d$var %in% Newjobs]            <- "New jobs"
  d$gruppe[d$var %in% Measures]           <- "Supportive measures"
  
  return(d)
}

d <- FillGroup(d)
d$order <- 1:nrow(d)

GetLabels <- function(d){
  offset <- c("   ")
  
  d$var[d$var %in% Enddate] <- paste(offset,c("2025","2030",
                                              "2040","2100"))
  
  d$var[d$var %in% Costs] <- paste(offset,c("EUR 0", "EUR 6", "EUR 12",
                                            "EUR 18"))
  
  d$var[d$var %in% Lostjobs] <- paste(offset,c("-5,000","-10,000","-15,000",
                                               "-20,000"))
  
  d$var[d$var %in% Newjobs] <- paste(offset,c("5,000","10,000","15,000",
                                              "20,000"))
  
  d$var[d$var %in% Measures] <- paste(offset,c("Expansion of renewables",
                                               "Support for new businesses",
                                               "Investment in modern infrastructure",
                                               "Research & development",
                                               "Training & early retirement"))
  return(d)
}

d <- GetLabels(d)

# bring in sublabels           
d$order <- 1:nrow(d)
dd <- data.frame(var= c("End date:",
                        " ",
                        "Annual costs:",
                        "  ",
                        "Lost jobs:",
                        "   ",
                        "New jobs:",
                        "    ",
                        "Supportive measures:"
),order=c(.5,4.1,4.5,8.1,8.5,12.1,12.5,16.1,16.5),
pe=1,se=1,upper=1,lower=1,gruppe=NA)
d <- rbind(d,dd)
d <-d[order(d$order),]
d$var <- factor(d$var,levels=unique(d$var)[length(d$var):1])

yylab  <- c("Phase-out Support (binary indicator)")

p = ggplot(d,aes(y=pe,x=var))#,colour=gruppe))
p = p + coord_flip(ylim = c(-.3, .3))  
p = p + geom_hline(yintercept = 0,size=.5,colour="black",linetype="dotted") 
p = p +geom_pointrange(aes(ymin=lower,ymax=upper,width=.4),position="dodge",size=.6)
p = p + scale_y_continuous(name=yylab,breaks=round(seq(-.4,.4,.2),1),labels=c("-.4","-.2","0",".2",".4"))
p = p + scale_colour_discrete("Attribute:") + scale_x_discrete(name="")
print(p)

theme_bw1 <- function(base_size = 13, base_family = "") {
  theme_grey(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.text.x =       element_text(size = base_size, colour = "black",  hjust = .5 , vjust=1),
      axis.text.y =       element_text(size = base_size , colour = "black", hjust = 0 , vjust=.5 ), # changes position of X axis text
      axis.ticks =        element_line(colour = "grey50"),
      axis.title.y =      element_text(size = base_size,angle=90,vjust=.01,hjust=.1),
      legend.position = "none"
    )
}

dev.off()
pdf(paste(ffilename,".pdf",sep=""),width=10,height=10) # for paper
p = p  + theme_bw1()
print(p)
dev.off()


